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1.0 INTRODUCTION 


The Interplanetary Program to Optimize Simulated Trajectories (IPOST) is 
intended to support many analysis phases, from early interplanetary feasibility 
studies through spacecraft development and operations. The IPOST output 
provides information for sizing and understanding mission impacts related to 
propulsion, guidance, communications, sensor/actuators, payload, and other 
dynamic and geometric environments. 

Much of the overall architecture for IPOST has been derived from the Program to 
Optimize Simulated Trajectories (POST) (Reference 1-1). Indeed certain POST 
parameters and capabilities have been incorporated into IPOST to aid in POST- 
IPOST user compatibility. IPOST has extended trajectory capabilities to target 
planets and other celestial bodies with intermediate and velocity correction 
maneuvers. IPOST capabilities and limitations are summarized in Table 1-1. 
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FEATURE 

CAPABILITY 

Optimization method 

Explicit (Master/subproblems), Implicit (collocation) 

Optimization algorithm 

NPSOL 

Optimization 

parameter* 

AV magnitude, mass, time, . . . 

Maximum controls 

25 (Master), 45 (subproblems), 1700 (collocation) 

Control parameters* 

Values of event criteria, AV, arrival conditions, thrust, 

Maximum targets 

25 (Master), 45 (subproblems), 1700 (collocation) 

Target parameters* 

Time, position, velocity, orbital conditions, . . . 

Targeting method 

NPSOL. Newton-Raphson, special Onestep 

Sensitivity matrix 

Finite differencing, analytic for special interplanetary 
targeting 

Maximum events 

100 

Event criteria* 

Time, distance, speed, closest approach, . . . 

Event activities 

Info, impulsive AV, launch, orbit insertion, mass 
jettison 

Maximum maneuvers/ 
subproblems 

15 


Conic, Onestep, Multiconic, Encke, Cowell, implicit 

Planetary bodies 

Sun, nine planets, Earth's moon, any user-defined 
bodies 

Ephemeris 

Analytic, precision (JPL) 

Trajectory 

perturbations 

Central body, perturbing bodies, radiation pressure, 
J2, aerodynamics, thrust 

Input/Output frames 

Ecliptic or planet equator, Mean 1950 or Mean 2000 

! * User selectable 1 


Table 1-1. IPOST Features/Capabilities 


IPOST, along with members of its family, such as POST and IPREP, can analyze 
and support almost every activity associated with space exploration. 

IPOST is event driven. That is, the user defines a sequence of events which are 
executed in the simulation process. The events can be triggered by different 
criteria, such as absolute or relative time, distance from a body, or propellant 
consumption. At the event times, various activities can be initiated or 
terminated, such as employing a different thrust steering law, changing 
trajectory propagators or propagation step size, performing an impulsive delta 
velocity maneuver or jettisoning a probe or stage. 

The time period between two contiguous events is called a phase. Trajectory 
propagation takes place in each phase. Five types of propagators are available 
(listed in order of increasing accuracy and decreasing computational speed): 
Conic, Onestep, Multiconic, Encke, Cowell. Propagator selection depends upon 
user needs, such as simple fast simulations for parametric feasibility analysis, or 
precision detailed trajectories to support subsystem design. 
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IPOST can run a single trajectory simulation or it can run multiple simulations. 
For multiple simulations, one can run a parametric scan and/or an optimization 
mode. The search mode will vary one parameter, such as planetary arrival time, 
over a specified interval and increment size, and perform a simulation (or 
optimization) for each search parameter value. 

The optimization mode will optimize a user cost/objective function, such as 
maximum mass that can be placed in a desired orbit, subject to user-specified 
constraints. The constraint variables, such as periapsis altitude or orbital 
inclination, are called dependent variables or target parameters. The parameters 
which are free to vary, such as maneuver delta velocity (AV), are called 
independent variables or control parameters. As part of, or instead of, 
optimization, trajectory targeting can be performed. In this case, there is no cost 
function and the IPOST problem reduces to finding a set of control parameter 
values that meet specified target parameter conditions. 

Generalized targeting and optimization uses the Stanford NPSOL algorithm. For 
certain types of problems, a trajectory decomposition method is available. There 
is a master optimization process which requires that the trajectory be divided into 
legs or sub-problems. Each subproblem is an optimization problem in itself,, 
containing controls, constraints and an (optional) objective function. A special 
application of decomposition is the Interplanetary Targeting and Optimization 
Option (ITOO). This technique uses analytical partials generated during nominal 
trajectory propagation to determine minimum AV (or mass) trajectories, usually 
for gravity assist (swingby) missions. 
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In addition to the classic method of explicit optimization, there exists an option to 
perform implicit optimization using the collocation method. In this case, each 
phase is divided into independent segments which are allowed to vary subject to 
intersegment continuity and the equations of motion. Optimization using 
collocation is less sensitive to faulty initial guesses, but requires much greater CP 
time than explicit optimization to achieve the same level of accuracy. 

IPOST input is via three namelists: $TOP, $TRAJ and $TAB. $TOP contains a 
description of the targeting and optimization problem. It must be input first. 
$TRAJ contains data that describes each mission event/phase. It must follow 
$TOP , and there must be one $TRAJ for each event. $TAB is used to input 
tabular data such as thrust vs. time or drag coefficient vs. mach number and 
angle of attack. Input and output units are metric. 
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2.0 COORDINATE SYSTEMS 

There are many types of coordinate systems used in mission analysis 
applications. IPOST provides a number of systems to allow the user a maximum 
amount of analysis insight and flexibility. 


2.1 INERTIAL ECLIPTIC SYSTEM 

This system is used during heliocentric (sun-centered) interplanetary flight, 
although it is sometimes used as a fixed reference frame for the entire mission. 


Center of system 

Primary plane 
X - axis 


the sun, planetary bodies, or 
planetary satellites 
Earth ecliptic 

Earth vernal equinox direction 



Figure 2-1. The inertial ecliptic coordinate system. 

(The center can also be at any of the planetary bodies and satellites) 
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2J2 INERTIAL P T.ANETARY EQUATOR SYSTEM 

This system is used during flight near a planetary body. 


Center - planetary body 

Primary Plane - planetary equator 

X - axis - rotation of planetary vernal equinox 

direction through right ascension 
and declination angles of the 
planetary pole vector. 



eq 


Figure 2-2. The planet equatorial coordinate system. 
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2.3 UVW SYSTEM 


This system is used for vehicles whose longitudinal axis or whose thrust axis is 
along the velocity vector. 


Center 

Primary Plane 
Primary axis (X) 


S/C center of mass 

orbital plane of S/C 

cross product of S/C velocity U = V x 

W direction and S/C angular 

momentum direction 



V 


— » 
v 



W 


-4 -4 

r x v 

-4 -4 

I r x v I 


U = Vx W 


Figure 2 - 3. The UVW Coordinate System. 
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2.4 RTN SYSTEM 

This system is used for vehicles whose longitudinal axis is along the local 
vertical, such as gravity gradient stabilized s/c. 


Center 

Primary Plane 
X - axis 


S/C center of mass 
orbital plane of S/C 
radius vector direction 



R 



N 


r x v 


I r x v I 


T = N xR 


Figure 2-4. The RTN Coordinate System. 
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2.5 THE ROPY FRAME COORDINATE SYSTEM 


The body frame is used in conjunction with other reference frames to orient the 
vehicle in celestial space, and to identify locations and orientations of vehicle 
components, such as the primary thrust vector and antenna boresights. 


Center 
x b axis 

z b axis 

y b axis 

P 

R 

Y 


- S/C center of mass 

- from center of mass through nose of S/C, 

- from center of mass through bottom of S/C, 
orthogonal to x 

A A 

- x b x z b 


pitch 

- roll 

- yaw 



Figure 2 - 5. The body coordinate system. 
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2.6 THE B-PLANE 


The B-plane coordinate frame is used for hyperbolic approach to, and departure 
from, a celestial body. It is often the most numerically stable system for 
describing planetary approach/departure conditions. 

HYPERBOLIC PATH 
OF SPACECRAFT 



0 = Orientation of B Relative to T 


S = Parallel to Incoming Asymptote 

T = Parallel to Reference Plane (Ecliptic Unless Otherwise Specified) 

— > — ^ — > 

R = S x T 


Figure 2-6. B-Plane Coordinate System 
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2.7 


The cone-clock system is used for determining the orientation of vehicle sensors 
and actuators. One application is to transform IPREP and QTOP thrust data into 



Figure 2-7. Cone/Clock Coordinate System 
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It can be seen from Figure 2 - 7 that to convert this vector to the heliocentric 
ecliptic (h/e) system, two rotations must be made. First, a rotation about xcc by 
the declination of the sun centered heliocentric ecliptic s/c vector, and then a 
rotation about the new z to line up the x and y axes. From this new vector in the 
h/e system, pitch and yaw can be calculated for the s/c orientation. This is input 
to IPOST. Since the cone clock system is not an inertial system, even if the cone 
and clock angles were fixed, the h/e acceleration vector would be continuously 
changing. Therefore, IPREP will calculate, for the IPOST user, cubic polynomial 
coefficients for yaw and pitch, as well as the acceleration magnitude, which can 
be input to IPOST as well. These polynomials are time dependent and can be used 
as controls in IPOST to "fine tune" trajectories due to modelling differences when 
higher orders of accuracy are desired. 

The orientation of the body system relative to one of the other systems is given 
through rotations of the following Euler angles: roll, yaw, and pitch. These three 
angles are defined by time dependent quadratics where 


roll 

= <t> = <t>! 

+ <)> 2 t 

+ <t> 3 t 2 

.[2-1] 

yaw 

II 

< 

II 

h- 1 

+ y 2 t 

+ y 3 t 2 

,[2-2] 

pitch 

II 

CD 

II 

CD 

+ 0 2 t 

+ 0 3 t 2 

,[2-3] 


The rotation matrices for body system to system for which roll, yaw, and pitch are 
defined as: 

"cos0 cosy cos0 cos<f> siny + sin0 sin<j> cos0 sin<J> siny + sin0 cos4> 

- siny cos(|) cosy cos(|) siny 

sin0 cosy sin0 cos<j) siny - cos0 sin<f> sin0 sin<{> siny + cos0 cosc|> 
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UVW system to planet equatorial system: 


U x 

V x 

W 

X 

IT 


W v 

y 

y 

y 

U z 

V z 

w 

z 


RTN system to planet equatorial system: 






L Z 


n 3 

N v 

n; 


z J 


Since each of these matrices consists of orthogonal rotations, the inverse rotation 
is merely the transpose of the matrix. 
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3.0 CONTC PROPAGATORS 


IPOST uses two classical propagators for two-body motion: GOODYR and 
LAMBRT. Both propagators can handle elliptical (including multiple 
revolutions) and hyperbolic orbits. 

GOODYR computes a final cartesian state given an initial state and a time of 
flight. The algori thm is based on Goodyear’s method (Reference 3-1). The two- 
body equations of motion are solved using a generalized eccentric anomaly and 
generalized Kepler's equation. The final state is computed using f and g series. 
The state transition matrix, and its inverse, are similarly computed. 

LAMBRT computes the initial and final velocities given the initial and final 
positions and a time of flight. This classical Lambert problem is solved using a 
method described in Reference 3-2. A generalized time of flight is iterated upon 
using an independent parameter which also depends on known orbit quantities. 



4.0 ONESTEP (1STEP) PROPAGATOR 


ONESTEP is a special adaptation of a generalized Multiconic method. It is used 
in JPL's PLATO, and described in References 4-1 and 4-2. The 1STEP 
computational sequence is as follows (see also Figure 4-1). 

(1) Start with an initial state at t^ and a desired propagation final time tf. 

There are other specifications needed, such as identifying the primary and 
secondary bodies and the sphere of influence time (At). 

(2) Estimate tp, the closest approach time at the secondary body, by Lagrangian 

interpolation using 3 radii (at initial, final, and midpoint times) based on 
primary body conic motion. 

(3) Compute the new corrected closest approach time tp' by conic propagation. 

(4) If tp and tp' are within an acceptable tolerance, proceed to Step 5, 
otherwise, replace tp with tp' and go back to Step 2. 



Figure 4-1. ONESTEP flyby process 
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(5) Propagate the initial state ti to t p using a two body conic (IPOST uses 
GOODYR) centered about the primary body. This state is called the 
(incoming) pseudostate. The propagation is from point 1 to point 2 in figure 
A-l. 

(6) Transform the pseudostate to the secondary body frame. 

(7) Propagate the state from tp to t p - At (backward in time) assuming the 
velocity at t p is constant. This is the rectilinear segment from point 2 to 
point 3. 

(8) Propagate the state at t p - At to t p + At using a two body conic (GOODYR) 
centered about the secondary body. 

(9) Propagate the state from t p + At backwards to t p assuming the velocity at t p 
+ At is constant. 

(10) Transform this (outgoing pseudostate) state to the primary body frame. t 

(11) Propagate the pseudostate from t p to tf using a two body conic (GOODYR) 
centered at the primary body. This is the final output state. The final 
propagation is illustrated as the segment from point 5 to point 6 in Figure 4- 
1 . 

There are also tests and modifications performed if the final time is before t p , or if 
the initial time is after t p , or if the trajectory is totally outside of the secondary 
body sphere of influence. 

The targeting and optimization process usually requires a state transition matrix. 
When 1STEP is used as a propagator, the state transition matrix can be computed 
analytically. This is done by chaining a sequence of separate transition matrices, 
all of which are computed analytically. 


$fi = aX/dXi = $65 ^54 ^43 ^32 °21 , [4-1] 


where $65 » $43 > $21 are con ^ c transition matrices generated as part of the 

GOODYR conic propagation in Steps 11, 8, and 3, respectively. The constant 
velocity transition matrices are 


^54 = '°32 


fl lAt\ 
K 1 ) 


, [4-2] 
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The Multiconic propagator is similar to ONESTEP (Section 4) except that it takes 
multiple steps to traverse a given time interval. This is done to avoid the 
restrictive assumptions of ONESTEP, particularly having to always fly to closest 
approach of the secondary body, and to accommodate additional perturbations. 
The algorithm solves the three body problem following References 5-1, 5-2, 5-3. It 
is non-directed, meaning that it does not presume a direction toward or away 
from either the primary or secondary bodies. It is presumed that the secondary 
body is in orbit about the primary body. 

The algorithm is straightforward. We are interested in propagating a S/C state 

[ Rj, Rj ] with respect to a primary body at time tj to another time, t F = tj + h in 
the presence of a secondary body. The solution is 

Rp = R ipf + r jgp * - hrj + (i (pp - Pj - h Pj) + 2 3ph ’ ^ ^ 


% ~ ^IPF + r ISF " r I + ^ (pF'Pl) + ha p ,[5-2] 

where Rjpp = position vector of S/C with respect to primary body mapped with a 
2 body conic from tj to tp. 

r ISF = P os iti° n vector of S/C with respect to secondary body mapped with 
a 2 body conic from tj to tp. 

r j = position vector of S/C with respect to secondary body at tj. 

Pj = position vector of the secondary body with respect to the primary 

body at tj. 


H = Hg/(p p +p s ) 

h = the delta time, tp - tj. 

a p = the average perturbing acceleration over At. 

The perturbing acceleration is the sum of the radiation pressure and the 
gravitational perturbations from bodies other than the primary and secondary 
bodies and other forces. Section 8 describes the perturbing acceleration models. 

There is an optional variable stepsize algorithm, 
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f. 


h = 


2 A 


S 0 R P r smin 


m 


^SMIN J 


where S Q is a user input scale factor, 


, [ 5 - 3 ] 


R t 


= distance of S/C from primary body, 


R SMIN = distance of S/C from closest secondary or perturbing body, 


^SMIN = gravitational constant of closest secondary or perturbing body 
The state transition matrix, if needed, is given by 

d* (tp, tj) = d>2B (tp, tj) + O 2B fp * o I J ’ 


Where 


<J> 2g is the two-body transition matrix from tj to t F about the primary body, 

g 

d> 2 b is the two-body transition matrix from tj to t F about the secondary body. 
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6.0 ENCKE PROPAGATOR 


For more precision trajectory propagation, the Encke method has been adapted as 
described in Reference 6-1. Instead of integrating the sum of all accelerations, as 
in the Cowell method, this method integrates the difference between the true orbit 
and an unperturbed conic orbit. Therefore, only the perturbations are integrated. 
Since the perturbations change much more slowly than the actual state, a larger 
step size can be utilized, making the method 3 to 10 times faster than the Cowell 
method, with a reduction in roundoff error. The reference conic orbit is the orbit 
which would result if all perturbing accelerations were removed at a particular 
time. 


The reference orbit is used for calculating the difference in total accelerations, 
until this difference becomes too large. When this occurs, a process called 
rectification occurs so that integration may continue without loss in accuracy. 
This means that a new epoch is chosen, and a new reference orbit is calculated 
which coincides with the true orbital path at that time. 

The objective of the Encke method is to find an analytic expression for the 
difference in the acceleration vectors of the true and the reference orbits. First, let 
denote the state vector of the S/C on the true orbit, and let t* denote the state of 
the reference orbit at the same time. The equations of motion for these two states 
are then 


and 


d 2 ? 

dt 2 


Ji 


— > 
r 



d 2 

dt 2 


a 

p 3 


— > 
p 


= o 


[6-1] 


[ 6 - 2 ] 


where "a is the vector sum of all perturbing accelerations. The difference vector 
of the true to the reference orbit is then 

5? = t - ? , [6-3] 
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and the second derivative with respect to time is given as 


or = r 


p 


[6-4] 


For initial conditions, recall that at the epoch, to , the state of the true and 
reference orbits are coincident, thus 


~r (t 0 ) = p* (t 0 ) and ~r (t Q ) = 'p > ' (t Q ) . 


[6-5] 


By substituting equations [6-1] and [6-2] into equation [6-4], we arrive at the 
analytic expression 


8 r = a p + 


p3 P 'r 3 r 


[ 6 - 6 ] 


which can be rewritten as, 


o — — > h- 

5 r = a p + p3 


[fi £\ 

— > £ — > 
r - o r 

Ll r J 

- 


[6-7] 


The above equation presents some numerical difficulties. The reason for using 
this method instead of the Cowell method was to obtain more accuracy, through 
less roundoff error. The expression 

( 1 - p 3 / r 3 ), with its very nearly equal terms, is difficult to evaluate numerically, 
requiring many more digits of accuracy. 

In reference 6-1 the problem was treated as follows. 

Let 


P 3 

-flq) = 1 - = 1 - (1 + q) 


3/2 


[ 6 - 8 ] 


where 


q=- 


(5 ~r - 2 r* ) » 5 7* 


[6-9] 


and 


3 + 3q + q 2 

fiq) = q i + (l + q) 3/2 


[ 6 - 10 ] 


Then, equation (6) becomes 


5 

5 r = 


-^3 [ tfq ) ^ + 5 "?] 


— ^ 
+ a P 


[ 6 - 11 ] 
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which removes the numerical problem, and is the equation integrated in the 
Encke method. 

The Encke method reduces the number of integration steps needed on a given 

interval since 8 ~r changes much more slowly than r . If the acceleration due to 

the perturbations approaches that of the central body term, or if 8 r / r does not 
remain small, the advantages of the method diminish rapidly. As the 
perturbations become quite large, then it is possible the reference parameters 
should be changed since the perturbations are becoming primary, or the Cowell 

method should be employed. In the case that 8 r / r has become large, the 
reference orbit more than likely needs to be rectified. 
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6.1 KNCKF, STATE TRANSITION MATRIX FORMULATION 

As an option in the Encke method, the state transition matrix, <J> , can also be 
integrated along with the state. 


From [6-1] we have 


ajfo) 

dlt (t 0 ) 


[ 6 - 12 ] 


where j^(t) is the state at time t, and x*(to) is the state at the initial time, to- 
When t = to, the matrix is a 6x6 identity matrix. The derivative of is given as 


where 


<D' = F<I> 

F _ dx^'(t) 
9 x*(t) 


[6-13] 

[6-14] 


If ~r is the radius vector, ~r is the velocity vector, and ~r is the acceleration 
vector, then 



Therefore we have the 6x6 matrix 


where 



I 

C 


<J) = 3x3 null matrix, 

I = 3x3 identity matrix, 


B = 


d~r 


and 


[6-16] 


[6-17] 


C=^ . [6-18] 

dr' 

For each set of perturbation equations, the matrices B and C are calculated, and 
used to create the matrix F. Once F is calculated for all perturbations, <&' can then 
be integrated along with the state. 

Perturbations which are accounted for in the transition matrix computation 
include perturbing bodies, J2, and atmospheric drag. 

6.2 ENCKE VARIABLE TIME STEP FORMULATION 

The Encke method also has the option of a variable time step which is dependent 
on the size of the perturbations. Reference 6-2 suggests that the time step should 
be dependent on the B partition of the F matrix derived in the preceding section. 
The time step for input to the numerical integrator (RUNGE4) is defined as 

h = ( I Bi 1 2 + | B2 1 2 + I B3 1 2 ).25 * TOL / SPDY [6-19] 

where 

Bx is the first vector column of B, 

B2 is the second vector column of B, 

B3 is the third vector column of B, 

TOL is the error tolerance level, and 
SPDY is the variable for seconds per day. 
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7.0 COWETX PROPAGATOR 


The Cowell method is the most straightforward of all the numerical trajectory 
propagation methods which include perturbations. The application of the method 
is simply to write the equations of motion of the studied object, including all 
perturbations, and then to integrate them step-by-step numerically. The 
equations of motion for a spacecraft would be 

~r" + (|i / r$) ~r = a^ [7-1] 

where a p is the vector sum of all the perturbative accelerations. This equation is 
broken down, for numerical integration, into a set of first order vector component 
equations. 


X, 

II 

X 

-> ' 
V X 

— ) 

= a px 

- (p / r3) x 

[7-2] 

II 

<< 

-» ' 
v y 

= a*py 

- (p/r3)y 

[7-3] 

I - > 

Z = V z 

-> ■ 
v z 

II 

N 

- (p/r3) z 

[7-4] 


Advantages of the Cowell method come from its simplicity of formulation and 
implementation. The advantages are equally weighted, though, by a set of 
disadvantages. When the acceleration forces are large (such as motion near a 
large attracting body), decreasingly smaller time steps must be taken which 
severely affect the speed at which the method operates. Also, when such small 
steps are taken, there can be a large accumulation of roundoff error. Roundoff 
error will also begin to accumulate in interplanetary flight if the step size taken is 
quite small. This gives support to the use of the Encke and Multiconic methods 
for interplanetary flight and for low perturbation trajectories. When the 
perturbations become quite large though, the Cowell method must be utilized. 
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8.0 PERTURBING ACCELERATION 


Perturbing accelerations include radiation pressure, low and high thrust 
propulsion, aerodynamics, gravitational nonsphericity of the central body and 
perturbations resulting from multiple attracting bodies. 
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8.1 RADIATION PRESSURE 


Acceleration due to solar radiation pressure is modeled by 


— ^ 
a 


SqCr 

m r^ 





) 

r , 


[8-1] 


where 

So 

Cr 


m 

-» 

r 



Solar flux constant. 

Coefficient of reflectivity, 

Mass of S/C, 

Heliocentric radius unit vector, 

Effective area used to calculate radiation press 

(orthogonal to body axes, transformed into the i 
frame). 


8.2 PROPULSION 


There are three types of propulsion systems that are modeled: low thrust, 
generalized, and blowdown. Low thrust propulsion accelerates ions that are 
supplied through electrical power systems, nuclear or solar. Generalized systems 
are traditional high thrust propulsion systems used in a variety of rocket 
applications, from launch vehicles to S/C attitude control. Blowdown systems are 
a relatively economical S/C propulsion system. 

All three propulsion systems provide both acceleration and propellant mass flow 
information to the propagators. The mass flow is integrated simultaneously with 
the translational equations of motion, using the mass flow equations that are 
applicable for each propulsion system. 

The thrust direction is assumed to be fixed with respect to the body. The thrust 
direction is given as a unit vector in the body frame. 

Acceleration due to low thrust electric propulsion is modeled as 


— ^ 
a 


2 u r[ P 
" m gl sp U ’ 


and the mass flow rate is given as 

m 1 a* 1 

m 

g lsp 

where 

u 

= Throttle level (normalized), 


= Thruster efficiency. 

g 

= 9.806 m/s, 
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[ 8 - 2 ] 



*cl 


I S p = Thruster specific impulse, 

= Thrust direction, 

= Electric power available to thruster. 


The available power for the thruster system is 


D D f, Cl C2 C3 ^ -P (t - to) p 

P = PO |^1 + r 2 r 3 + jA J e L ' P HK’ 


where 


PO 

Ci,C 2 ,c 3 

PL 

to 

PHK 


= Initial power, watts, 

= Solar array efficiency factors (=0 for NEP), 

= Reactor decay factor (=0 for SEP), 

= Initial time, 

= Housekeeping power for space vehicle, watts. 


Two finite thrust models will be employed. The first is a generalized system 
which can be a regulated pressure system in which a separate gas tank provides 
constant pressure to the propellant tank, or a solid rocket. The generalized finite 
thrust is modeled as 

T = uT vac - A E™ > 


T -> 
a = — u 
m 


where 


I or , = Thruster specific impulse, 

bp 

m = Mass flowrate of propellant to thruster(s), 
u = Throttle level (normalized), 

ii = Thrust direction, 
m = Mass of the S/C, 

T, ro . = Vacuum thrust of the engine(s), Newtons, 

V db 

A £ = Exit area, m 2 , 

P(h) = Atmospheric pressure at altitude, h. 
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A blowdown system will also be modeled. A blowdown system uses one 
pressurized tank partially filled with propellant. The system is modeled as 


m = u m 


max’ 


a = 


T 0 u 11 


m 


1 + 


u m (t - to) 
p f V 0 Ul 


Y ’ 


[8-5] 


T = ml a 


I =T/(m g), 
ht 
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where 
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m 

m 

u 

— > 
u 

m 

Pf 

to 

Vo 

Ul 
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max 


= Initial thrust level, Newtons, 

= Mass flowrate of propellant to thruster(s), 
= Maximum mass flowrate, 

= Throttle level (normalized), 

= Thrust direction, 

= Mass of the S/C, 

= Propellant density, 

= Initial time, 

= Initial tank volume, 

= Ullage ratio, 

= Ratio of specific heats. 


8.3 GRAVITATIONAL PERTURBATIONS 

Disturbing body accelerations are modeled in the Multiconic, Cowell and Encke 
propagators as 


— > 
a 


= £ - -> i 


[ 8 - 6 ] 


i=l I r • - R I 


where 


= Gravitational constant of the i^ 1 body, 

= Position vector of the i^ 1 body (from ephemeris), 

= Position vector of the space vehicle, 

= Number of perturbing bodies. 

Gravitational potential for nonspherical J harmonics is modeled as 


Pi 
— > 
r i 

-4 

R 


n 
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where 


(I = gravitational constant for the central body, 

L = planetocentric latitude, 

r e = equatorial radius. 


30 


Therefore, gravitational accelerations due to J 2 in the planet centered frame are 
given as 


2 

[8-7] 


[ 8 - 8 ] 


[8-9] 

8.4 AERODYNAMICS 

Aerodynamic accelerations are composed of drag and lift, with no side force. The 
drag and lift coefficients can range from constants to bivariant functions of mach 
number and angle of attack. 

The atmospheric model used in aerodynamic calculations can be either a simple 
exponential density or pressure/temperature versus altitude. If the exponential 
atmosphere model is selected, certain aerodynamic effects will not be included, 
such as thrust back pressure and mach number variations. 

Acceleration due to aerodynamic drag is modeled as 

t -i C D & p v a V a , [8-10] 
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8x x 


_8<j) _ 
: 8z - 


p z r 
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where 

C D 

A 

m 

P 

v a 


= drag coefficient, 

= cross sectional area of S/C perpendicular to the 
direction of motion 

= S/C mass, 

= atmospheric density at the altitude of the S/C, 

= I ~r & I = speed of S/C relative to the 
rotating atmosphere, 
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where 



r 


W + 6'yA 
y' - 9'x 

v z ' J 


6c 


k z j 


inertial velocity in planet 
equatorial system, 
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where 

Cl = lift coefficient, 

n ' = unit vector in the lift direction. 

As mentioned earlier, two types of atmospheric models are available for 
calculation of the density at a given altitude. The first is a simple exponential 
model where the density is 

(-l/ah)(h - h 0 ) 
p = Poe > 

where 

p 0 = base density, 

a^ = scale height, 

h = current altitude, 

ho = base altitude. 

Currently stored in p 0 , a^ , and h 0 are the values for the Earth's atmosphere. 

Also available is a more complex tabular atmospheric model which calculates, for 
a specific altitude, atmospheric density, atmospheric pressure, molecular weight, 
mach number, and molecular temperature. Stored tabular values are altitude, 
molecular weight, molecular temperature, and pressure. These values are 
interpolated to give the correct output values for density, pressure, weight, mach 
number and temperature. As with the exponential atmosphere, the stored values 
currently available are for the Earth. 
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9.0 ANAT.YTTC EPHEMERIS 


The analytic ephemeris in IPOST is taken from MAPSEP (Reference 9-1). Nine 
planets with respect to the sun are represented. Each planet is characterized by 
six conic elements (a,e,i,o>,£l,M) in a heliocentric mean ecliptic 1950.0 frame of 
reference. Each conic element is described by a third order polynomial in time. 
For example, Mercury's inclination in radians is 

i = .01222... + (-1.041...E-4)*T + (1.745...E-8)*T 2 + (0)*T3 

where time (T) is measured in Julian centuries past 2000.0 . 

When the ephemeris routine (EPHEM) is called to compute the planet state 
(position and velocity), the conic elements are evaluated via their time 
polynomials. The conic elements are then transformed to cartesian space, and 
the state vector is output. 

Analytic planetary satellite ephemeris is treated similarly. For example, the 
Earth’s moon is described by six conic elements with respect to Earth, each of 
which has four polynomial coefficients describing the respective element in time. 
At any specified time, the lunar cartesian state is computed, and added to the 
Earth cartesian state to provide the moon’s state with respect to the primary body. 
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10.0 P RECISION EPHEMERIS 


Extremely accurate planetary and lunar positions and velocities are derived by 
evaluating a file of Chebyshev polynomials calculated/tailored for a given mission 
application. This is accomplished using multiple-day-arc Chebyshev polynomials 
whose coefficients are derived by a pre-processor using ephemeris data supplied 
by the Jet Propulsion Laboratory (JPL), Pasadena, California. These data files 
contain positions and velocities referred to the rectangular equatorial coordinate 
system of the mean equator and equinox of 2000 for the major solar bodies: 
Mercury, Venus, Earth- Lunar barycenter. Mars, Jupiter, Saturn, Uranus, 
Neptune, Pluto, and Earth's moon plus the nutation rates in longitude and 
obliquity. This data was generated at JPL by weighted, least-squares estimation of 
appropriate orbital models using source positions obtained on the basis of current 
planetary theories. 

Using the method of Reference 10-1, a solar/lunar/planetary (S/L/P) file of 
Chebyshev polynomial coefficients are generated for the time period(s) of interest 
for each mission application. For trades performed for given mission(s), these 
files would be saved for later access to avoid recomputing the coefficients. In this 
manner, Chebyshev polynomial representations will retain the accuracy of the 
original JPL data while increasing computing efficiency by eliminating the need 
to interpolate on the JPL ephemeris data tapes. The Chebyshev file accessed 
during computations contains polynomial coefficients for each component of 
position and velocity and for each element of the matrix that transforms from the 
mean equator and equinox of date to the true of date coordinate system as required 
for IPOST calculations. Also included in the file would be the coefficients for the 
equation of equinoxes, AH, used to correct mean Greenwich sidereal time. 

In conjunction with the ephemerides pre-processor, a file-read utility is required 
to access the Chebyshev polynomials and to calculate the ephemerides for IPOST 
or to compute them in the coordinate system(s) of choice for output. 
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11.0 NPSOL 


NPSOL is a sophisticated optimization package developed by Systems 
Optimization Laboratory of Stanford University (References 11-1 and 11-2). It 

minimizes a smooth cost function F(~y ,~u) where ~u are the control (independent) 
parameters subject to upper and lower bound constraints on both the independent 

variables and on dependent functions, ~y{u). These ~y can be linear and/or 
nonlinear functions, some of which may be target (dependent) variables. 

NPSOL is a stand alone package which interfaces with IPOST by accepting current 
trajectory values (the cost function F, the cost gradient ~g = dF/d~a , the cost Hessian 
H = d 2 F/d~u 2 , controls *u , constraints y , and the sensitivities C = d y Id u ) and 
returning the next iteration control parameters u which minimize the cost 
function F("u ). This process is iterated to convergence, eventually outputing the 
optimum and the values for F and y . 

The NPSOL algorithm (Reference 8) uses major and minor iterations to solve the 
targeting and optimization problem. In the major iterations NPSOL seeks a 

significant decrease in the merit function along a direction of search p . NPSOL 
defines the merit function as: 

F("?) = I [ li * (ci (^) - si)] + |l [ri * (dfi?)- si) 2 ] 
i i 

where F("u ) is the objective function 

li is the vector of lagrange multipliers 
ci is the i-th constraint function 

si is the set of slack variables used to handle inactive constraints and 
ri is the penalty vector for constraint violations. 

Therefore, during NPSOL operation, the problem becomes better targeted and 
more optimized in a simultaneous fashion. 
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In the minor iteration NPSOL seeks the search direction P for the major iteration 
by minimizing the quadratic programming problem: 

Q = ‘gT ? + !?TH?, 

subject to a set of constraints, where ~g is the gradient of F( u ), and H is the quasi- 
Newton approximation to the Hessian of the merit function. 

Both PGA and NPSOL are gradient based algorithms. They require derivatives of 
the objective function and of the constraints with respect to the control variables. 
The former vector of derivatives is called the objective gradient and the latter 
matrix of derivatives is called the Jacobian, or sensitivity, matrix. In addition, 
NPSOL constructs a second derivative matrix or Hessian. 

One of the keys to optimization success is a well conditioned Jacobian matrix. 
Usually, this is formed by finite differencing. In IPOST, the finite difference 
control perturbations can be input by the user or computed by NPSOL. If NPSOL 
computes the perturbation size for each control parameter, the appropriate 
perturbation is that control value which produces the most accurate partial 
derivative for the control-constraint (or objective) combination, which is just above 
the n um erical noise threshold. It can use either forward or central differencing 
techniques. Quite often about a third of the rim is spent computing an accurate 
Jacobian because many simulation passes, or function evaluations, are needed to 
produce the finite difference perturbation. NPSOL also has the option of 
"verifying" an input perturbation size by constructing and comparing a finite 
differenced Jacobian. 
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12.0 TRAJECTORY DECOMPOSITION 


A complete, end-to-end mission analysis for complex interplanetary missions 
with one or more intermediate flybys is not practical in one computational step. 
The complexities arising from gravity-assist and velocity change maneuvers 
present a formidable challenge to solve only one trajectory through the 
deterministic boundary value problem. Inclusion of all perturbations required for 
optimization and their subsequent cascade of multiple trajectory calculations 
would expand the problem to an unmanageable scope. Decomposition of the 
mission into tractable sub-problems has long been a viable alternative (Reference 
12 - 1 , 12 - 2 ). 
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12.1 TCXPTJCTT OPTIMIZATION - MASTER/SUBPROBLEMS 


Decomposition is the resolution of the mission into sub-problems which are 
readily solved, followed by an iterative master level solution of approximate sub- 
problem solutions. The problem is described as a series of events and trajectory 
phases, specified by user input. The events can reflect physical happenings 
(starting time, time of periapsis, impulsive velocity change) or be informational 
events (time to record the state vector or a sphere-of-influence). Each trajectory 
phase starts and ends on an event. Complex missions are decomposed by 
partitioning the trajectory into subproblems of one or more phases. Sub-problems 
usually mirror a natural grouping of trajectory phases - such as a planet to planet 
leg. 

Each subproblem is defined by a set of control (independent) parameters and a set 
of target (constraint or dependent) parameters. The controls are varied such that 
target conditions are met. This requires first generating a target sensitivity 
matrix, or Jacobian, ususally by finite differencing. In IPOST there are two ways 
to adjust the controls using the Jacobian matrix. One option is to employ a 
Newton- Raphson technique to "invert" the Jacobian. This is done with a L-U 
decomposition algorithm which does not require a full rank system, that is, the 
nimber of controls and constraints do not have to be equal. The second option 
employs NPSOL. NPSOL can be used to optimize a subproblem unique objective 
function at the same time meeting constraints. An example of subproblem 
controls would be midcourse AV, and subproblem constraints would be planetary 
encounter conditions at a subsequent flyby. The master level problem supplies the 
starting independent variable values to each subproblem leg and receives in turn 
the targeted values. 

The master problem, which represents the complete mission design, is the 
summation of the individual sub-problems. It controls optimization and specifies 
the individual sub-problem targets. The sub-problem solutions are relatively 
small in dimension. The master problem iterates through the sub-problem 
trajectory legs to optimize an overall target strategy in conjunction with the 
NPSOL optimization program. NPSOL evaluates the total trajectory in terms of 
target error and modifies the independent variables for the next master level pass. 

12.2 IMPLICIT OPTIMIZATION - COLLOCATION 

A more subtle form of decomposition is applied in the technique of collocation 
(Ref. 12-3, 12-4). This reformulation of the optimization problem requires setting 
up the trajectory as an implicit simulation where the trajectory is divided into 
many "independent" segments. 

Implicit simulation computes the vehicle state between two junction points, or 
nodes. This trajectory segment can be a portion of a simulation phase or the 
entire phase. At nodes between events, the vehicle state must be the same for 
each side, that is, just prior to and just after the node. At nodes that correspond to 
events, the states on each side may be different. This discontinuity would occur if 
instantaneous event activities were specified, such as mass jettison or an 
impulsive AV maneuver. 


39 


The vehicle state is represented between nodes as Hermite, or cubic, polynomials. 
Each state component corresponds to a single Hermite polynomial. To compute 
the polynomial coefficients requires position, velocity, and acceleration (from 
evaluating the equations of motion) at each node. Figure 12-1 describes the 
methodology and compares it with explicit propagation. 

In collocation, the states on each side of each node, that is, the pre- and post-node 
states, are allowed to vary. These free states become additional control variables. 

Additional constraints must also be introduced because the states have other 
imposed conditions. For example, the state on each side of an event node is 
related by the event activity. Thus, an additional constraint becomes the node 
connectivity between pre and post states, particularly due to any state 
discontinuity arising from a mission activity, such as an impulsive AV. 

Another set of constraints is introduced because the state on each side of a node 
which is internal to a simulation phase, or between events, must be the same. In 
some collocation formulations, the zero state differences for these internal nodes 
are introduced as constraints. However in IPOST, the state on each side of an 
internal node is explicitly set equal to the state on the other side. This has the 
effect of reducing the number of controls and constraints such that only a single 
state at each internal node become control variables. 

One other IPOST refinement eliminates the state differences of the first event 
node by explicitly setting the post-event state equal to the pre-event state plus any 
state discontinuity associated with activities in the first event. In theory, all 
double states can be reduced to single states by this procedure, but this remains a 
future modification. 

There are other constraints introduced by collocation. These constraints involve 
forcing the Hermite polynomials to meet the equations of motion at the mid-point 
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(tc ) of each segment. The constraints, called defects, ((?) are formulated as: 

• • 

cf = X (X, tc) - X pjp _ o 

x = [ x , *y , z , x , y , z , m ], the seven element vehicle state 

derivatives 

where X (X , tc ) = equations of motion evaluated at t c using Hermite polynomials 
(x) evaluated at 

V 

X HP = analytic derivations of Hermite polynomials, e.g., 
if y = bg +bit + b2t 2 + b 3 t^ 

then yjjp = b i + 2 b 2 + 3b3t 2 


The net result of collocation is that for any reasonable simulation/optimization 
problem there are hundreds more control and constraint parameters than in an 
explicit simulation/optimization(see Figure 12 - 2). In the IPOST formulation, an 
equal number of controls and constraints are added via collocation. For example, 
if a classical problem has 10 controls and 6 constraints, and the simulation has 5 
phases with 4 segments per phase, then the total number of controls and 
constraints are : 

Number of controls = 10 + nc * 7 = 192 

Number of constraints = 6 + nc * 7 = 188 
where nc = (number of segments per phase + 1) (number of phases) + 1 
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controls: x = state 

■* ... 
constraints: E = event state (dis)continuity 

-4 

d = state derivative defect 




-4 

E 

I 

event 

node 


| ^ | 

1 ' t 

internal internal 

node node 


"d 


% 


t 


1 


4 4 

X .X, 


■4 

E 


4 


internal event 

node node 


Figure 12-2. Controls and Constraints Introduced by Collocation 

Although the dimensionality of the optimization problem has been greatly 
increased, there are several compensating features to collocation which 
ultimately expand its region of convergence. The first benefit is that each 
segment is relatively independent of the others because of the free states at each 
node. Thus, a change early in the mission is not amplified into extremely large, 
and unpredictable, changes later in the mission, as in explicit 
simulation/optimization. 

Another collocation benefit is that most derivatives are analytic, and therefore 
rapidly, and accurately, computed. The few derivatives that must be done by finite 
differencing are of the "local" variety. That is, there is very little time mapping, at 
most over a single segment, which makes the finite differenced partials relatively 
stable and well-behaved. 

Finally, the total Jacobian, although quite large, is a sparsely populated, banded, 
matrix. There are many techniques for efficient generation and manipulation of 
sparse matrices. In IPOST, the zero elements of the Jacobian are identified and 
pre-set, thus eliminating the wasteful application of finite differencing by NPSOL. 
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13.0 INTERPLANETARY TARGETING AND 
OPTIMIZATION OPTION (ITOO) 


The optimization objective is defined as a function, F( u ), which is to be minimized 

over some time interval, t 0 to tf , where u is an m x 1 set of control parameters 
(independent variables). 

The targets are defined as a set of n x 1 constraint parameters (dependent 
variables) 

~y(u) = ~y (n , x(u), c ), [13-1] 

where ~x (u ) = x* (u), x (t 0 ) are state variables (position, velocity, and perhaps 

— ^ 

mass) and c is the array of constraint values on y . 

The general targeting and optimization technique is a recursive application of "u 

new = old + a"u , where a"u depends on Ay* (a measure of how far away a target 

is from its desired value or constraint boundary) and the gradients dF/d u and 

d~y/d~a . The Au solution may also require past information on a'u, A"y , and/or the 
gradient matrices, plus numerous user-supplied weighting factors. 

The problem solution starts with an initial guess of u 0 and x 0 = x (t 0 ). The 
equations of motion are integrated from to to tf (or some user defined stopping 

condition that defines tf). The A~y and F functions are evaluated. The cost and 
constraint gradients are computed. A control correction, a^j, is computed and the 
new u is formed. The new conditions, and x* 0 , are propagated to form a new 
trajectory. This process is repeated until F is minimized (i.e., until dF/d"u = 0). 

Assume that 

F = I ,wiyi 2 , [13-2] 

i=l 

where L is the number of variables of interest, wi is a weighting factor, and yi is 
an unconstrained dependent variable and 

f = ?(V,V). [13-3] 

From a truncated Taylor series expansion, 
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[13-4] 
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This provides a direction of A u and an initial estimate of its magnitude. This 
technique has special applications that will be shown later. It presumes that y is 

"small" and that *y and "u are linearly related. These conditions are quite often 
true in interplanetary situations. 

The local Jacobian at time t is defined as dz (t)/9x (t), with ~z representing various 
trajectory parameters, such as B-plane parameters, and x the state parameters 

(position and velocity). ~z can be control or target variables. The local Jacobian is 
computed by numerical differencing. Although explicit analytical partials exist 
for many interplanetary parameters, future growth dictates the need for the more 
generalized numerical differencing method. 

There is a broad class of interplanetary problems which lends itself to be solved 
using the methods of ITOO. Figure 13-1 illustrates tlfe situation. The initial and 
final times are specified. There are N bodies that are flown past. An impulsive 
trajectory correction maneuver is performed between each body. The final body 
often contains terminal conditions, such as planetary injection or small body 
rendezvous. Intermediate planets usually have constrained flyby conditions. 



Figure 13-1. General Interplanetary Mission 

Note that this geometry has applications beyond the transplanetary phases and 
includes satellite tours, such as the Galilean satellite tours, cycler missions 
that visit the same two bodies, such as Earth-Mars cyclers, and small body tours, 
such as multiple asteroid missions. 

The targeting/optimization problem decomposes into two stages. The first stage is 
a simple targeting problem. The impulsive AV of each maneuver is computed 
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such that the following body's encounter conditions are met. There are three 
components of A V and three target variables (user selected). 

The second stage contains the optimization. The function to be optimized is 
identical to that discussed earlier. Here the dependent variables, y , are the 
maneuver A 'Vs. The independent variables, or control parameters, "u , are the 
maneuver times, Ti , and the planetary encounter conditions, Pj , which are given 
as 
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[13-6] 


Here we note special requirements affecting the initial and final conditions. In 
transplanetary applications, there is typically a launch planet, P 0 , whose 
launch/escape conditions can constitute an additional set of control parameters. 
Also, there can be one or more terminal maneuvers at the terminal body, such as 
rendezvous or orbit capture, which must be included in the total cost function. 

A reasonable upper bound for the number of bodies appears to be 15. This is 
consistent with all missions currently being projected for this century. The 
n um ber of control variables is 4n+3 and the number of dependent variables is 
3n+3. Hence, the maximum number of control and dependent parameters should 
be 63 and 48, respectively. In reality, IPOST limits control and dependent 
parameters to 45 each. 


There are no constraints on the a Vs because they are part of the cost function to be 
minimized. A few words on cost function weighting are in order. 

If all the wis are unity, the cost function is simply the sum of squares of the AV 

magnitudes. A more useful cost function might be the sum of the total AV 

46 


rn 


magnitudes. This can be done by setting wj =1/1 AV. I . An even more useful 

function would be the sum of the total propellant expended (which is equivalent to 
maximizing delivered mass). This can be done by using the rocket equation, 


wi = 


mi 


I AVI 2 



-IaV l/Ci 


[13-7] 


where mj is the S/C mass prior to AV , and Cj is the exhaust velocity and is given 

by 


Ci - So^sp-* 


[13-8] 


There are, or can be, constraints imposed on the control variables. Typically, P Q is 
constrained by launch vehicle conditions, such as C 3 or declination. Maneuver 
times are usually restricted from occurring near a planet so as not to interfere 
with science data collection. Planetary encounter conditions are also constrained, 
such as not allowing closest approach beneath the planet's surface (impact). 

Thus, we have a cost function dependent solely on unconstrained dependent 
variables, subject to constrained independent/control variables. 


Whatever the optimization method, a necessary part is the gradients, 8F/(3u ) and 

d"y/d"vL Because of the cost function formulation, the two gradients are related as 
follows 



n +1 — > T 
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[13-9] 


We need only to compute dA\^/d u and we have both cost and target gradients. 

The computations are even more efficient if we make use of analytic partials. 
Onestep, Multiconic and Encke propagators can produce analytic state transition 
matrices as part of their trajectory propagation process. 
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The state transition matrix is composed of four partitions, 
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^j-l.i-l , ^i j-l, Ojj,... are all output from the Onestep and Multiconic propagation 
process. We also note that 

^i, i-1 — ^i, j-l ^j-1, i-1- [13-11] 

The local Jacobian can also be easily computed, and made available as 
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Subsequent chaining can generate other needed partials by 
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dx . dx . 

J J 

The first step of the targeting/optimization problem is targeting. A simple 
iterative scheme can be employed where 

d~V. . 
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new A V . = old AV . + . 
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[13-14] 


where 
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and APj is the error in the encounter conditions (desired minus actual) from the 
last iteration. The partials, d~V :/d~P should be updated for each iteration. 

Once targeting of each leg has been completed, the second stage is minimization 

of the cost function. Again, the chain rule can be applied. The results of 3A V/9u 
are summarized below. 
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For the sensitivities of AV. to maneuver times, Tk, where i=l,...,n+l and 

k=l,...,n, 

1) for k > i, 
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= 0 


[13-21] 


which means that all planets beyond the i^ 1 target planet have no effect on dAV . ; 
2) for k = i, 
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[13-22] 


[13-23] 


4) for k £ i-2, 


dA^V. 

W 


= -SyR2( i ) 


_ n s RR1 (L) 

L=k+1 


A^ 


For illustration, the full expression for k = i-2 is written 


dA*V. 

i 

fav. 

1 

d~R. 

l 

di t. 

i 

^,1 

3Ti-2 = ‘ 

— ^ 



_ ~ > 

dR. 

^ 1 

3R . 1 
L l-l 

3v i-i 

dR. , 
1-1J 


— ^ ^ ^ 

\ 

— > — > 

dV. dV. dV. , 

1 1 1-1 



dR. , dR. , dV. „ 
l-l l-l i-2 

+ — > — > 




dR. , dV. , dR. , 
l-l l-l 1-1J 

J 


dR- 2 av._ 2 dR._ 2 _ 



[13-24] 


[13-25] 
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Similarly, we can compute the sensitivities of AV. to planetary encounter 
conditions where i=l,...,n+l and k=0,l,...,n, 

1) for k > i, 


dA~V. 

A 


= 0 ; 


2) fork = i, 


3a"v 

d~P. 

1 


1 • SvpoO) ; 


3) fork = i-l, 


dA~V . 

A 


3R . 3V. 

1 1 


SvRO^) " 

. i i 3v.. av.,i 

i-l L i-l 1-1J 


SyPoG'D > 


4) for k = i-2, 


-v ^ ^ 

3AV. 3r. 1 

1 =S V R 2 «)^T :J - Svpo(i-2) 


d~P. 


i-2 


d'v. 


i-2 


5) fork^i-3, 


dA~V . 


a?, 


1 = S VR2^) 


n s RR1 (L) 

L=k+1 


alt 


^ Svpo(k) 


av 


[13-26] 


[13-27] 


[13-28] 


[13-29] 


[13-30] 
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Again, for illustration, the full expression for k — i-3 is written 


dAV . 


dR: 3R. 


BR. d~R. dV. , 
11 l-l 

— + ^ _ 7T _ 

)R i-i 3v i.i 3R i.i 


— ) — ) "4 

av. av. av. , 

1 1_ w 

— > + — > — > 

3Rj.l 3 v i-i 3R i-l 


^ -v ^ ^ 

8R. ] ^ 3R; .[ dVj.2 

3R i-2 37 i-2 3R i-2 


3R i-2 37 i-3 
3 ^-3 3R i-3 


Wc have now computed all of the partials needed to construct 3F/3 u and 3 y /3 u . 

Optimization may cause one or more of the planetary legs to become untargeted. 
Thus, the targeting and optimization stages may have to be repeated until 
satisfactory convergence and optimization is achieved. The methodology of two- 
stage interplanetary targeting and optimization is consistent with the IPOST 
trajectory decomposition capability. 
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14.0 SPECIAL ONESTEP TARGETING 

TARG1S is a restricted application of 1STEP which is a special adaptation of a 
generalized Multiconic method. This specialized targeting scheme is suitable for 
a restricted class of interplanetary problems. (See references 4-1 and 4-2 for 

restrictions.) An impulsive change in velocity ( A V ) is computed at a specified 
maneuver time (T 0 ) to achieve target body relative conditions (targeting) at 
periapsis. Target conditions which work well are radius of closest approach 
(RCA), B -plane orbit angle (0), and time of closest approach (T p ). It basically 
solves a two point boundary value problem (TPBV) between the pseudostate 
position and desired target conditions (see Figure 14-1). In theory, this is more 
robust than solving the TPBV between initial velocity and target conditions 
because the nonlinear effects of the primary body propagation from initial state to 
pseudostate are minimized. 



primary body 


Figure 14-1. Special ONESTEP targeting 
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Define E as the target conditions at the desired time Tp. The target controls 
currently implemented are listed in the User's Manual. A change in the 
pseudostate position, R*, will affect E in two ways: directly and indirectly (through 


V*). 


5E = 


f 


BE 


BX 


BX BX . 

y p SOI 

V 


— ^ ^ 

ax . ax, 

SOI 

ax„ 


Be 


Bx 


B~b L 


Bx Bx . 

p SOI 


^ ^ ^ -s ^ A 

ax . ax* 3 v* av 

SOI _ 0 

^ ^ ^ -v ^ 

ax* av*av o 3 r*^ 


V_ 




J 


direct 


indirect 


We can compute J = 


-» -4 

Be/Bx 


the Jacobian by finite differencing. 


Furthermore, 




= <£>p S is the state transition matrix from Tsoi to Tp from GOODYR. 


ax 


BX 


so - = <J> S * is the constant velocity transition matrix = 


'I -IAT" 
O I 


BX, 
3 R* 


and 


ax, 

a\^ 


av, 

av 


= lower right 3X3 portion of O* 0 from GOODYR. 


a'V 
c 

B~R, 


lower left 3X3 portion of <!>„< from GOODYR. 


3 V* 3V o ^ /T ^ ^ N (l - AAT^j 

If we define A = — zr > then = (J) (^sp) I A I 

av 3R„. 3R * v J 


— ^ ^ 

av 3r 

o 


* v AV * 

The targeting procedure is as follows. 

— >T ( — ^ ^ 

(1) Start with initial primary body centered state X Q = I R Q V q 
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and desired secondary body centered end conditions E(J. 


(2) Propagate X to T p with respect to the primary body to obtain 
xT = p* V*) and 

(3) Change pseudostate to secondary body frame of reference, and propagate 
back to T S oi = Tp - AT where AT = secondary body sphere of influence, 
with constant velocity, then propagate to Tp with respect to the 

secondary body to obtain X and Op S . 


— > — > 

(4) Evaluate the actual target conditions, E , from X . 

a p 

(5) If the target error AE = Ed - Ea meets desired tolerances, then exit 
targeting (go to Step 11). Otherwise, continue targeting (go to Step 
6 ). 

(6) Compute the Jacobian, J = 9 e /d X . 


(7) Form ~E/d~R *, as expressed above. 


(8) Compute S = 1 d~E/d R 


(9) The new R * = old R * + S AE . 


(10) Solve Lambert's problem from R to R * with flight time of Tp - To . 

This provides a new initial velocity V . Go to Step 2 and continue 
the targeting process. 

(11) The total impulsive AV is the last Vq - original Vq . 
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15.0 SPECIAL FUNCTIONS 


15.1 LAUNCH MANEUVERS 


The launch maneuver sub-event in IPOST applies a maneuver from a parking 
orbit to a hyperbolic escape orbit. Computed outputs from this special subroutine 

are the new S/C mass (after escape burn), the AV applied to complete the required 
maneuver, the bum duration, and the state of the S/C at periapsis of the 
hyperbolic escape orbit. 


a = - = semi-major axis of the hyperbolic orbit where Voo is the 

v 2 

oo 

magnitude of the outgoing asymptote vector of the hyperbolic orbit. Voo is either 
input or is calculated from 


= VC3 


A 


V„ 


cos 0 

COS 

a 

cos 8 

sin 

a 


sin 8 


where 

5 = desired declination, 
a = desired right ascension, 
C3 = desired escape energy. 
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Periapsis radius Rp = hp + R s 
where 

hp = circular park orbit altitude, 
R s = surface radius. 



A 

T 


V x z 

oo 

IV x z I 

oo 


A A A 

R = T x V 
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If i is input, then 


A 

A 

A 

B = 

COS 0 T + 

sin 0 R 

A 

A 

A 

R = 

-cos 6 V 

+ sin <J) B 

P 

T oo 


A 

A A 


N = 

B x V 



oo 


A 

A A 


V = 

N X R 


p 

P 



b V M 


II 

R 



P 
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Outputs: 


NOTE: 


ENDNOTE. 


where 



escape state at periapsis 


If F/(m g) < .02, then a low thrust spiral out is assumed and 

av = V% + v £ 


m F 
At = 


- (AV/(g*Isp)) 

= m 0 e 

mo AV 
F 


Isp = specific impulse 
F = average thrust during maneuver 
g = reference gravitational acceleration 
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15.2 IMPULSIVE MANEUVER 


The impulse maneuver sub-event in IPOST applies a maneuver, such as a 
midcourse correction, to the S/C state. Computed outputs from this subroutine 
are the mass of the S/C after the maneuver, the bum duration, and the in and out 
of plane AV angles. These are calculated as follows: 

AV = I a"v I 

- (AV/(g*Isp)) 

nip = mo 6 

where 

mp = final S/C mass after the maneuver, 
m 0 = initial S/C mass before the maneuver, 
g = surface gravitational acceleration of the earth, 

Isp = specific impulse. 


Am = mp - mo 


where 


At 


mo AV 
F 


= bum duration, 


F = average thrust during maneuver. 


A 


Orbit normal N = 


R x V 


I l?x ~V I 


A 

AV = 


— > 
AV 

AV I 


A A 


P = cos‘1 (AV • N) gives 0< P ' < k 


**• P = 


n 


P gives 


n 


< p < 


7t 

2 


where 


P = out of plane AV angle 
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A A A A 

T = N x (AV x N) 


A A A , 

a' = directed angle from V to T about N gives 0 £ a <2 k 

( If a '< 7 i then a = a' ^ gives - n < a < n 

I If a ' > n then a = a ' - 2 tc J where a = in plane AV angle 
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15.3 ORBITAL INSERTION 


The orbital insertion maneuver sub-event in IPOST takes a S/C on a hyperbolic 
flyby trajectory and performs a maneuver to insert it into an orbit about a planet. 
Calculated outputs are total AV required for the maneuver, final mass after the 
maneuver, total bum duration, state of the S/C at periapsis of the new orbit, and 
the time from the current state to the final orbit periapsis. These are calculated 
as follows: 


Rph = a (1 - e) 


where 

a = semi-major axis for the hyperbolic orbit, 
e = eccentricity for the hyperbolic orbit, 

Rph = radius of periapsis for the hyperbolic orbit. 

Rp = hp + Rs 

where 

hp = desired periapsis altitude, 

Rs = surface radius of the planet, 

Rp = radius of periapsis of the desired orbit. 

Ra = ha + Rs 

where 

ha = desired apoapsis altitude, 

Ra = radius of apoapsis of the desired orbit. 

Rph + Ra 
a l =_ 2 


Ra + Rp 

a 2 = 2 

where 

a^ = semi-major axis of intermediate orbit after first AV, 
= semi-major axis of the desired orbit. 


Vph = 




2a - Rph 
L a Rph J 


velocity at hyperbolic periapsis 
before the first maneuver, 
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= velocity at hyperbolic periapsis 
before the first maneuver, 



If iD = desired inclination < 0, and Ra < Rph, then 



where 

i = inclination of the hyperbolic orbit. 
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Otherwise, 


AV X = Vph -V ai 

m 0 AVi 
At l “ F 


= m 0 e - 


(AVi } 
g*Isp 


where 

F = average thrust during the maneuver, 
= mass after the first maneuver. 




2 n 

2a^ - R a 

a l Ra 

L .. 

2^2 • Ra 


velocity at apoapsis of transfer 
before the second maneuver, 


velocity at apoapsis of desired 
after the second maneuver. 


If iD = desired inclination < 0, and R a < Rph > 

av 2 = ^v p 2 i + V a 2 j -2VpiV ai 


then 


cos (iD - i) 


where 

i = inclination of the hyperbolic orbit. 

Otherwise, 


AV 2 =V Pl -Va 2 


m l aV 2 


At 2 - p 


m 0 = mie 


f av 2 

g*Isp 


\ 

) 


where 
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F = average thrust during the maneuver, 

= mass after the first maneuver. 

The transfer time, TF, from the hyperbolic periapsis to the desired orbit periapsis 
is given by 



The total AV is therefore given as 
AV = AVi + AV2 . 


NOTE: If the value [ F / m 0 g ] < .02 then the maneuver is a low thrust spiral in, 
and therefore 


A V = V (p/Rp) - p/a) 
TF = At = m 0 AV / F 


mf = moe 


(AV/g*Isp) 


a 2 = 

If mean anomaly, M, is greater than zero, then M = -M, 
t p = -M/n + TF 

where 

n = Vp T- a?. 
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Compute the final orbit state at periapsis, x , 


— > 


Convert orbital elements with M = 0, to position, r , and 
velocity, v . 


where 


Therefore 


R p = (Rp / Rph ) r 
V p = (V p / V p h ) ~v 


Vp = V (2a2-Rp) /a2Rp) 


— ^ 
X 


ft ^ 

p 
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15.4 SPACECRAFT MASS 

The N-stage mass properties model from POST has been adapted for IPOST. The 
mass of the spacecraft (m s / c ) after an event (positive side) is calculated from the 
mass before the event (minus side) as 

mj c = m g/c - (mjett + m P r °p) ' Am 


where 


mjett = 



mprop = 
= 

Am = 


mass to be jettisoned, 

a factor to specify the fraction of propellant mass to be 
dumped, 

current (remaining) propellant mass 
propellant use (scale) factor 

change is propellant mass due to propulsive accelerations 
(Section 8.2) 


Finite thrust propulsion algorithms model low thrust, generalized thrust, and 
blowdown systems. These models return a mass flow rate which is integrated to 
obtain Dm using a simple Euler integration. Before each integration step the 
propellant (mp r0 p) is tested to determine if sufficient fuel exists for the next step. 

Impulsive maneuvers (launch, impulse, and orbital insertion) directly change 
mg/c (Sections 15.1, 15.2, 15.3, respectively). 
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